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We engineer a quantum bath that enables entropy and energy exchange with a one-dimensional 
Bose-Hubbard lattice with attractive on-site interactions. We implement this in an array of three 
superconducting transmon qubits coupled to a single cavity mode; the transmons represent lattice 
sites and their excitation quanta embody bosonic particles. Our cooling protocol preserves particle 
number—realizing a canonical ensemble— and also affords the efficient preparation of dark states 
which, due to symmetry, cannot be prepared via coherent drives on the cavity. Furthermore, by ap¬ 
plying continuous microwave radiation, we also realize autonomous feedback to indefinitely stabilize 
particular eigenstates of the array. 


Ordinarily, uncontrolled dissipation destroys quantum 
coherence, but it is now appreciated that an engineered 
quantum bath is a valuable resource for quantum compu¬ 
tation (state preparation p]-[4] and system reset(5[) and 
quantum error correction [Bj. A dynamical bath can in¬ 
duce cooling or heating, and is of great utility in optome¬ 
chanical ground state preparation 00, quantum state 
transfer between various EM modes Emu, amplifica¬ 
tion [12H - 114J . many-particle quantum simulation mm, 
and entanglement generation mm- 

Initially conceived and implemented in trapped ion 
systems [13 ED], dissipation engineering has been im¬ 
plemented in solid state systems with two recent exper¬ 
iments on superconducting qubits: autonomously con¬ 
trolling the orientation on the Bloch sphere of a sin¬ 
gle qubit 0, and stabilizing a Bell-state in a two-qubit 
system [4j. In this Letter, we experimentally demon¬ 
strate that dissipation engineering can be used to control 
a novel, complex superconducting system embodying a 
much larger Hilbert space: the ten lowest-lying energy 
levels of a coupled, three-transmon array which realizes a 
one-dimensional, attractive Bose-Hubbard Hamiltonian. 
The techniques employed in this work define a path for 
cooling and stabilizing complex quantum systems to spe¬ 
cific target states. 

The Bose-Hubbard Hamiltonian [2T] is a prototypi¬ 
cal model used to describe a broad class of quantum 
matter. While the repulsive side of the Bose-Hubbard 
phase diagram has been extensively explored in ground¬ 
breaking experiments with ultracold atoms in optical lat¬ 
tices [22M24] . the attractive regime has thus far eluded 
emulation. Our realization of this model here, in per¬ 
haps its simplest incarnation, opens the door to exper¬ 
imental verification of as-yet unexplored predictions of 
attractive Bose-Hubbard dynamics: the existence of self¬ 
bound states [251127] and the possibility to create large- 
scale multipartite entanglement [25]. 

The cooling protocol we develop here, based on Ra¬ 
man scattering processes, facilitates entropy and energy 
exchange between the qubit array and its bath while pre- 




FIG. f: (a) A schematic (not to scale) of the experimental 
setup described in the text, (b) Spectroscopically measured 
eigenfrequencies of the one- and two-particle states of the ar¬ 
ray as a function of current through the external bias coil. 
For a given current, the flux through the two SQUIDs in the 
array differs by 2.5%; 17 mA roughly corresponds to a quar¬ 
ter of a flux quantum. Solid lines denote measured frequen¬ 
cies with fits to the ID Bose-Hubbard Hamiltonian shown as 
overlaid dashed lines. Red lines correspond to two-particle 
states; green lines are one-particle states. The inset shows 
raw data near the \Ei) frequency, from which the darkness 
of the |G) —> \Ei) transition discussed in the text becomes 
apparent. 

serving the total number of excitation quanta in the ar¬ 
ray. Essentially, this amounts to simulation within the 
canonical ensemble; a similar path to grand-canonical 
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FIG. 2: (a) and (b) Examples of spontaneous decays from 
two-particle states to the global ground state via the one- 
particle subspace. In (a), |Fs) decays rapidly and almost en¬ 
tirely to |i? 2 ), which then decays to |G), while in (b), | Fo) 
decays with roughly equal rates to \Ei) and I-E 2 ). Overlaid 
black lines are obtained from fitting the decay data for all nine 
excited states to a single rate-equation model, as described in 
the supplement [25| . (c) An illustration of the natural decay 
pathways of the system, from the two-particle subspace on 
the left, through the one-particle subspace in the middle, to 
the zero-particle state on the right. Each black arrow repre¬ 
sents a decay channel, with the opacity of the line indicating 
the rate of the transition. Darker lines indicate faster rates, 
as the legend shows. Also shown are representations of the 
eigenstate wavefunctions in the qubit basis. The black circle 
radius is proportional to the mean particle number. 


simulation via bath engineering (using photons) has re¬ 
cently been proposed theoretically by Hafezi et al 530]. 
In another related work, a cooling scheme similar to 
ours has been employed recently to measure the dynamic 
structure factor of a gas of cold atoms |3T. Addition¬ 
ally, a recent experiment with superconducting qubits 
has achieved simulation of the unitary Fermi-Hubbard 
dynamics via discrete gates [32] . 

Our system is comprised of an array of three 
capacitively-coupled transmon qubits [33], each coupled 
dispersively f34] to a waveguide cavity [35]. The two 
qubits on the ends of the array utilize a SQUID loop 
to allow tuning their frequencies via an external mag¬ 
netic field. Taking into account the full transmon spec¬ 
trum [33] . the system is described by the Hamiltonian 
H = H c av -f- H arra ,y “t“ Hint , with 


ffcav = (a) a + 1/2) (1) 


Harray = h V Ubfa + Sfbtybfo) 

+HJ J2 Q>]+ibj + b^-bj+i) 

j=l 

+bJi 3 (b\b 3 + 6361 ) 
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Hint = ft^gjjbjCL^ + b\g), (3) 

3 =1 

where and are creation operators for cavity pho¬ 
tons and excitation quanta of the j th qubit in the ar¬ 
ray, respectively. H cav is the Hamiltonian for the 7.116 
GHz waveguide cavity, which couples dispersively to each 
qubit in the transmon array with strength gj, as de¬ 
scribed by Hi n t. The array itself is described by U array : 
each transmon is a weakly anharmonic oscillator with 
|0) —> |1) transition frequency ujj and (negative) anhar- 
monicity ay. The three-qubit system is effectively an 
array of lattice sites on which particles—the transmon 
excitation quanta—can hop with nearest-neighbor tun¬ 
neling strength HJ , and next-nearest-neighbor tunneling 
hJ \3 (set primarily by the capactive coupling [36]), with 
J \3 < J. In this language the negative anharmonic- 
ity of the transmon gives rise to an attractive pairwise 
interaction between these particles, since distributing a 
pair of excitations among two identical transmons takes 
an energy a more than lumping the quanta together on 
the same qubit. It is this combination of tunneling and 
on-site interaction which realizes the canonical ID Bose- 
Hubbard Hamiltonian m- 

Since the Bose-Hubbard Hamiltonian conserves total 
particle number, eigenstates of the three-qubit array can 
be grouped into manifolds characterized by this quantum 
number. In our experiment, we work with the zero, one, 
and two-particle manifolds, comprising respectively one, 
three, and six states. We denote the zero-particle state 
by |G), the single-particle states by {| Ei), i £ [1,3]}, 
and the two-particle states by {|F}), j £ [1,6]}, with in¬ 
creasing subscript value indicating higher-energy states. 
Due to the dispersive interaction between qubits and the 
cavity, a reflected microwave signal near the cavity reso¬ 
nance frequency acquires a phase shift dependent on the 
state of the array. Amplification of the readout signal 
via a Joseplison parametric amplifier m and subsequent 
higher-temperature electronics allows for measurement of 
the signal phase and hence the array state. 

Using this dispersive readout, we first characterize the 
system by spectroscopically probing its energy levels as 
the edge qubit frequencies are tuned down via external 
magnetic flux. As shown in Fig. [l] our measurement of 
the system’s one- and two-particle energy states agrees 
well with predictions based on the attractive ID Bose- 
Hubbard Hamiltonian. The extracted parameters are: 
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FIG. 3: (a) An approximate representation of the eigenstates of the array-cavity system; states in the left-hand (right-hand) 
column contain zero (one) cavity photons. In this picture, our cooling process can be understood qualitatively as follows (taking 
\E 2 ) —> |£7i)) as an example: the cooling drive at frequency w c — A c induces a transition from IF 2 } |0) to |Fi) 11), where the 
second ket indicates the cavity photon number. The cavity state |1) decays via photon emission, leaving the system in the 
state |2?i) 10), as desired, (b) Example of cooling from \E 2 ) —> \Ei) at an incident power corresponding to 1.3 drive photons in 
the cavity, (c) The cooling rate versus drive frequency is Lorentzian, centered around lo c — (E 2 — E\), with linewidth roughly 
k. As the drive power is increased, the Lorentzian peak shifts due to a Stark shift of the relevant transition frequency. The 
inset shows that in the regime where the cooling rate is much lower than k, the rate scales linearly with incident power, (d) 
Example of cascaded cooling, where the system is cooled from F 3 ) to |Fi) via the intermediate state |F 2 ). (e) Connected graphs 
representing the measured cooling rates per photon in the linear regime for the one- (left) and two- (right) particle subspaces, 
with the width of the line indicating the rate of the corresponding transition. The dispersive shifts for the I-F 3 } and IF 4 ) states 
are almost identical, so they cannot be distinguished by our measurement. We thus do not measure cooling from IF 4 ) to |Fi). 


wi/2tt = 5.074 GHz, w 2 /2tt = 4.892 GHz, oj 3 /2tt = 5.165 
GHz, J/2 tt = 0.177 GHz and J 13 /2it = 0.026 GHz, all to 
within ± 0.003 GHz. The cavity-qubit coupling gi/2n = 
149 ± 7 MHz, g 2 /2-K = 264 ± 7 MHz, 53 /2tt = 155 ± 7 
MHz and qubit anharmonicities cri,3/27r = -214 ± 1 MHz, 
a 2 / 2 - 7 T = -240 ± 1 MHz, were calibrated independently. 
While at zero applied flux the qubits are relatively well- 
separated in frequency, as the edge qubits are tuned down 
towards the middle qubit, avoided crossings become ap¬ 
parent, showcasing the coupling between the qubits. Our 
system lies in the parameter regime where the compet¬ 
ing tunneling (J/2n ~ 180 MHz) and on-site interactions 
(a/2ir ~ —220 MHz) have nearly equal strength. 

At an external field corresponding to 10 mA in the 
bias coil, the individual qubit frequencies approximately 
coincide—the detunings between neighboring qubits, 45 
and 115 MHz, are lower than the tunneling rate J. 
At this bias we characterize—in the absence of engi¬ 
neered dissipation —the particle-loss dynamics resulting 
from coupling to the dissipative environment comprised 
of both the leaky cavity (n/2ir = 10MHz) and micro¬ 
scopic material imperfections. A coherent microwave 
pulse initializes the system to the desired | Ei) or | Fj) 


state, after which the resulting populations are measured 
as a function of time. As expected, these decay dynamics 
fit well to a model which only includes single-particle loss 
events: direct transitions from the fq) states to |G) are 
suppressed, as recently observed with a single transmon 
qubit [55] . Examples of such decays are shown in Fig. 
b, while Fig. H shows the full map of decay rates within 
these subspaces. 

A striking feature of the natural decay dynamics is the 
discrepancy between decay times of different states in 
the same manifold. For example, \E 3 ) and \E 2 ) live for 
~ 30/zs, while | E 3 ) decays much more quickly, in ~ 3/rs. 
This is due to the substantially different dipole transi¬ 
tion matrix elements that each single-particle state | Ei) 
exhibits with respect to the final state |G). A related con¬ 
sequence of these dipole moments is shown in the inset 
of Fig. [I] where at 10.71 mA, | E\) goes fully dark, i.e. 
becomes impossible to excite via a coherent microwave 
pulse. Under only Purcell decay- in the absence of ma¬ 
terial losses in the system—such a dark state would live 
indefinitely, making it attractive for shelving an excita¬ 
tion if it could be readily prepared. We will return to the 
preparation of these dark states as one application of our 
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FIG. 4: Steady state population at the different stages of 
the persistent stabilization scheme of the two particle ground 
state |Fi>. The top trace shows the thermal equilibrium pop¬ 
ulation, where 78% of the population is in |G). In each sub¬ 
sequent trace a further drive is added: first a coherent drive 
|G) —> £ 3 } (I), then cooling from I.E 3 } —> | Ei) (II), coherent 
drive from \E\) —► \Fi) (III), and finally a cooling drive from 
\F±) —> |Fi) (IV). At the end, 67% of the population is in the 
desired state, |Fi). [£ 4 ) contains the bulk, 13%, of the resid¬ 
ual population. Single-particle (|£i)) stabilization achieves a 
population of 80% (data not shown). 


bath engineering protocol, which we now discuss. 

By altering the quantum noise spectrum of the bath 
interacting with the Bose-Hubbard chain, we can mod¬ 
ify the system decay dynamics. In our cooling protocol, 
this alteration takes the form of an additional microwave 
drive incident on the cavity, red-detuned from the cav¬ 
ity resonance by an amount A c , as illustrated in Fig. 
for the \E 2 ) to |I?i) transition. This induces quantum 
photon shot noise whose spectral density peaks at fre¬ 
quency |A C | f35]. When HA C matches the energy dif¬ 
ference between the array’s initial state |i) and a lower 
eigenstate |/), the emission of a photon at the cavity’s 
resonance frequency mediates a so-called cooling transi¬ 
tion from \i) to I/). The energy gained from the cooling 
transition augments the incident photon energy to allow 
emission on resonance. Similarly, a blue-detuned drive 
induces heating transitions to array states of higher en¬ 
ergy. The rate for these processes will depend on both 
the Raman-transition matrix element between the initial 
and final array states and on the incident photon flux. 
The transition rate increases with photon flux up to a 
value of ~ ft, saturating there as the dissipative process 
requires the emission of a photon by the cavity. Since 
k is much larger than most of the natural decay rates, 
this technique is well suited to driving otherwise inacces¬ 
sible particle-number-conserving transitions in our sys¬ 
tem. More details on the cooling protocol can be found 
in [221. 


To characterize the cooling dynamics, we initialize the 
system into an | Ei) or \Fj) eigenstate and subsequently 
apply a cooling drive for a variable time and measure the 
state of the array. Cooling rates are extracted via a fit to 
a model similar to that used for the natural decays, with 
additional parameters to capture the induced intraman¬ 
ifold transitions. Because the cavity’s density of states 
exhibits a Lorentzian profile with width ft, so will the 
transition rate as a function of cooling drive frequency, 
as shown for the \E 2 ) to \E{) transition in Fig. [3]^. 

For incident powers which cool at a rate r,^ much 
lower than «, a Fermi Golden Rule calculation [22] shows 
that 


Ti-J./ OC Pin I Mif I 


(Wi - + A c ) 2 + (ft/2) 2 


(4) 


where Pi n is the incident cooling drive power and Mif 
the matrix element connecting the states |i), |/) of the 
cooling operator jj29]J which describes the effect of the 
dissipative bath; in other words Mif quantifies the cou¬ 
pling between the states |i) and |/) indirectly via the 
cross-Kerr terms that couple the qubits with cavity. The 
predicted linear scaling of the peak T with P- m is shown in 
the inset of Fig. [sjc for the \E 2 ) —>■ | E\) transition. |M^| 
provides a measure of the efficacy of each transition; we 
map out this value for each pair of eigenstates, showing 
the results in Fig. H- 

In most cases, applying a drive whose frequency is tar¬ 
geted to cool |i) to I/) has no effect on the decay dynam¬ 
ics of the other states, as most cooling drive frequen¬ 
cies are spaced apart by more than several ft. However, 
when multiple cooling frequencies are separated by less 
than ft, a single drive can give rise to a so-called cas¬ 
caded cooling effect, whereby multiple cooling transitions 
happen in sequence. In our system, for example, the 
I-F 3 } —► IF 2 ) and the IP 2 ) —> |Fi) transitions are sepa¬ 
rated by only 17 MHz, so a single tone can cause the 
system to cascade from F a ) to |Fi) via the intermediate 
state IF 2 ), as shown in Fig. § 1 . In the specific example 
shown, since the cooling matrix element between \Fz) 
and IF 2 ) is substantially lower than that of \F 2 ) and |Fi) 
{\M 21 \ ~ 5|M 32 |), we cooled with the drive frequency 
tuned to the I-F 3 ) —► \F 2 ) transition; this achieved ap¬ 
proximately equal cooling rates from |£ 3 ) to \F 2 ) and 
\F 2 ) to |Fi). Cascaded cooling sequences could be useful 
in larger many-qubit systems with manifolds containing 
several closely-spaced eigenstates (see [29] for details). 

The transitions |G) —> \Ei), |G) —> \E 2 ), and | Ei) —> 
|Fi) do not interact strongly with the electromagnetic 
environment of the cavity on account of the symmetry 
of the states; this decoupling is responsible for their rel¬ 
atively long lifetimes. Correspondingly, however, it is 
difficult to coherently initialize these states via direct 
transitions, but our cooling scheme affords their efficient 
preparation. To illustrate this, consider the |G) —► | E\) 
transition, which as shown in Fig. [T] is at its darkest at a 
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flux bias of 10.71 mA. At this bias point, we use the Ra¬ 
man cooling protocol to prepare | E\) indirectly via the 
I-E 3 ) —> | Ei) cooling transition (data not shown). Fur¬ 
ther, by combining coherent drives with Raman cooling, 
we stabilize \Ei) indefinitely against particle loss. As the 
first part of Fig. [4] illustrates, this is done by coherently 
driving the |G) to I.E 3 ) transition with a Rabi frequency 
of 7 MHz while applying a drive to cool the \E 3 ) —1 | E\) 
transition at a rate of 3 MHz. We next use \Ei) as a 
stepping stone to stabilize the two-particle ground state 
|Fi), as shown in the lower part of Fig. [4j To accom¬ 
plish this we add two additional drives, an extra coher¬ 
ent drive from | E\) —> \F^) with a Rabi frequency of 7 
MHz, and a cooling drive from \F 4 ) —y |Fi) with a rate 
of approximately 3 MHz. Observed fidelities, while in¬ 
line with a rate matrix calculation, are primarily limited 
by spurious thermal population of dark states, which can 
be reduced by additional cooling tones. This initializa¬ 
tion and maintenance of the array in the ground state 
of a specific particle-number manifold will be a valuable 
resource for a hardware simulator. 

In conclusion, we have realized a three-qubit trans- 
mon array and spectroscopically verified that it obeys 
an attractive ID Bose-Hubbard Hamiltonian up to the 
ten lowest-lying eigenstates, highlighting the use of cir¬ 
cuit QED in simulating otherwise challenging quantum 
systems. Our developed cooling and stabilization proto¬ 
cols, based on quantum bath engineering-a well-studied 
phenomenon in quantum optics-afford effective control 
over this solid state system. The capabilities demon¬ 
strated here-engineering decay dynamics and stabilizing 
particular eigenstates-slrow that dissipation engineering 
can be a valuable tool as superconducting circuits scale 
up in complexity to complement simulators based on cold 
atoms and trapped ions. 
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Supplemental Materials: Cooling and Autonomous Feedback in a Bose-Hubbard chain 

with Attractive Interactions 


THEORY OF COOLING IN A TRANSMON ARRAY COUPLED TO A CAVITY 

Here we describe the theory behind the bath-engineering protocol implemented in the experiment. We first consider 
the general case of an array of L transmon qubits, each capacitively coupled to its nearest neighbors. After showing 
how, in such an array, an engineered bath can induce transitions between array eigenstates which conserve excitation 
number, we specialize to the case of a three-site array (L = 3) used in our experiment. As we will show, theoretical 
predictions are in qualitative agreement with experimental observations, with the single-mode approximation of the 
3D microwave cavity likely contributing a large part of the discrepancy. 


General L-site array 

A general array of L capacitively-coupled transmon qubits (j = 1,... L) is described by the Hamiltonian 


L 

^array b ^ ' 

3 = 1 


ujjb^bj 




^ + hj 


L -1 

E 

j=i 


(bLrb 


'j + l u J + b ] b 3 + 1 )- 


(SI) 


The first sum describes each individual transmon as an anharmonic oscillator with creation/annihilation operators 
and bj. This is valid in the limit in which we work, where the Josephson energy of each qubit dominates the 
capacitive charging energy by roughly two orders of magnitude. Capacitive coupling, or dipole-dipole interaction, 


between the qubits gives rise to the second term in Eq. (SI), with J being the hopping amplitude for an excitation 


to jump from one qubit to the next. Here, for simplicity, we ignore couplings beyond nearest-neighbors, because the 
dipole interaction scales inversely with the cube of the qubit-qubit distance, making next-nearest-neighbor coupling 
amplitudes close to an order of magnitude lower than corresponding nearest-neighbor values. 


For a derivation of Eq. (SI) via circuit quantization, see the supplementary materials of either ;S36 ] or I S40j . These 
derivations specialize to the lowest two transmon levels, but considering higher levels is straightforward. In this work, 
we work in a parameter regime where the qubit-qubit coupling is much lower than the qubit frequencies; in such a 
regime the rotating-wave approximation is well-justified, and we use it here to obtain the Bose-Hubbard Hamiltonian. 
This Hamiltonian has U(l) invariance and therefore preserves the number of excitations. 

The array is coupled to a single-mode cavity H cav = fio; c a^a (with and a the creation/destruction operators) via 
an interaction term H lnt = 9 j(bjd^ + feja). The coupling gj between qubit j and the cavity depends on the 

product of the strength of the electromagnetic field of the cavity mode at the location of qubit j and the transition 
dipole moment of that qubit, and hence varies from qubit to qubit. For example, in a setup such as the one used in 
the main text, where the array lies roughly in the center of the cavity, the qubit-cavity coupling to the lowest cavity 
mode will be largest in the center of the array. 

Since the operators which create and destroy transmon excitations satisfy the bosonic commutation relation [bj, frj] = 
1, it is natural to view these excitations as bosons on a lattice, with each transmon embodying a lattice site. Each 
transmon has a negative anharmonicity ay < 0, representing an effective on-site interaction between the bosonic 
excitations. In the case where the transmon anharmonicity is uniform across the array, close to what we realize in 


the experiment, the array Hamiltonian in Eq. (SI) has the familiar form of a Bose-Hubbard model in the regime of 


attractive interactions (here a plays the role of the Bose-Hubbard parameter U). 

To complete the description of this circuit-QED system we add the drive term, representing an external classical 
drive applied to the cavity: Hdrive = he(t)(a) e~ luldt + ae lUdt ), where e(t) is the strength of the external drive and 
the drive frequency. The full Hamiltonian of the system qubit-array and cavity is the sum of all the aforementioned 
terms: 


H = H; 


array 


T ft c :av -f- Hint T ffflrive • 


(S2) 


This Hamiltonian describes all aspects of the array-cavity system important to our experiment, including the cooling 
via engineered dissipation, as we will show shortly. To see this, we cast the Hamiltonian in a simpler form via a few 
standard transformations. First, we move to the rotating frame of the drive by applying the unitary transformation: 
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U _ e *Kw d t(a t o+j] i bth 3 -)^ t rans f ormec j Hamiltonian H = UHW — iHUdtffl |S44j takes the form: 

L L -1 L 

H = + °<^bjb^bjbj'j + hJ y^(bt +1 bj + b]bj+i) + HA c ata + tC^g^bjO* + b]a) + fie(t)(a^ + a), (S3) 

i=i j=i j=i 

where we have introduced the detunings: A c = w c —and Aj = ojj — u>d- In essence, this rotation gives new effective 
qubit and cavity frequencies while removing the oscillating time-dependence on the cavity operators in the drive term. 

Ignoring for a moment both the drive and the anharmonic transmon terms in H , the array-cavity system is simply 
a set of coupled harmonic oscillators which exhibits its own set of normal modes. To find these modes, we write the 
quadratic terms of the Hamiltonian, H 0 , in matrix form: 

L L -1 L 

Ho = hy^ Ajb\bj + h.J E(^+i bj + b\bj + 1 ) + hA c a^a + fi^ g^bja) + b\a) (S4) 

j-i j= i j=i 

= Hv^Hqv, (S5) 

where we have introduced the L + 1-component vector v = (a, bi, 62 , • ■ • bt) consisting of the annihilation operators 
for each oscillator, and the (L + 1) x (L + 1 ) matrix "Ho representing their frequencies and couplings. Now the normal 
modes of the coupled system are found simply by diagonalizing the matrix "Ho- We thus find the matrix M whose 
columns are the eigenvectors of Ho] then M~ 1 HqM is diagonal. Calling N = M _1 , the corresponding change of basis 
W = iVv gives the eigenmodes of the system. Writing these modes explicitly in terms of the matrix elements of N 
and M will be useful for later purposes—calling the new basis vectors W = (A, Bi, B 2 , ■ .. I?l), we have: 


and their inverse transformations: 


A = N 00 a + ^ Noibi , 
1=1 

L 

Bj = Njoa + Njibi, 


1=1 


a = M 0 qA + M 0 iBi , 

;=i 

L 

bj = MjoA + y ( MjiBi. 


1=1 


(5 6 ) 

(57) 


(5 8 ) 

(59) 


The corresponding eigenvalues or normal-mode frequencies are denoted by AIn this new basis Hq is, as expected, 
a sum of uncoupled harmonic oscillators: H 0 = HXqA^A + h A jB^Bj. 

In the dispersive limit, where gj <C |w c — wj|, we expect the dressed mode A to largely consist of the cavity mode, 
with small contributions from each qubit. Equivalently, this means that M 0 i and N 0 i are each much smaller than 
M 0 o an( I Wo, which are both close to unity. Conversely, the Bj modes are mostly linear combinations of transmon 
excitations ( 61 , 62 , ■ • - 6 l) with a small component from the cavity. In symbols, |Mjo| <C \Mji\ and IW'ol ^ \Hji\ for 
j,l 0 (see for instance the specific example from the experiment, shown in Eq. (S50) and ( S51[ >). 

In this new basis the drive term becomes: e(t)(a^ + a) = e(t) ^Mqo(AI + A) + J2iL 1 M 0 i(Bj + 5/)^. Because of 


the mode mixing, the drive on the original cavity-mode a now excites all of the normal modes {A, {I?j}}. However, 
since M 0 i M 0 0 , we can neglect the Bj terms in this operator unless the drive is close to resonance with one of the 
qubit-like modes. In that case, the coefficients Moi determine how responsive the qubit-like mode is to a drive pulse, 
i.e., which states are dark and which are bright. We will return to this point later. 

Continuing onwards in solving the original Hamiltonian, we now put back the anharmonic transmon terms which 
we neglected when finding the linear eigenmodes. In terms of the new normal-mode operators, the anharmonic term 
is 



1 L L L L L 

2 E + E M H B \)( M j 0 rf + E MjpBl){Mj 0 A + Y, M jq B q )(M j0 A + ^ M JS B S \ S10) 

j —1 Z=1 p= 1 gr=l s=l 


Invoking the rotating wave approximation to neglect terms which do not conserve excitation number, the remaining 
terms can be grouped into three categories: 
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• self-Kerr corrections to the qubit-like operators: J2i pqs =i MpqsB]B^BqBg, 
where we have defined the tensor: 

L 

f-^lpqs — E OljMjiMjpMjqMja, (Sll) 

3 = 1 

• a self-Kerr correction for the cavity: HqA' A' AA, 
where we have defined the constant: 

L 

n 0 = E/ a j M jO' (Si2) 

3 =1 


• and a cross-Kerr term that couples together the A and Bj operators: 4A AJ2i p Vi p B}B p , 
where we have defined the tensor: 


L 

1 lip = E a jMjoMjiMj P 
3 =1 


(S13) 


Now we have the full Hamiltonian in the new basis W: 

H w = hXoA^A + A4 + h ^Z X 3 B j B 3+2 E HipqsBjBpBqB s + 2 hA^ A ^ ] r li P BjB p 

j =1 Zp<j's=l Ip 

L 

+he(t)M 00 {A f + H) + fie(f) ^ M ol (B} + B t ). (S14) 

;-t 

Given our assumption |Mjo| -C |Mjj| for l,j ^ 0 we have: no <C r/i p <C Hi pqs - In other words, in this new dressed-state 
basis W we have: [qubit self-Kerr] [cross-Kerr] [cavity self-Kerr]. For simplicity we henceforth neglect the cavity 
self-Kerr term n qA^A^AA. 

Next, we displace the new cavity-like operator A and write it as a classical part and a small quantum correction: 

A = A(t) + D, (S15) 

by applying the unitary transformation: Ujj = exp [A(t)A^—A*(t)A\. Under this transformation Ur>AU' D = A—A(t) — 

D and the Hamiltonian Hw transforms as: 

H w ,d = U d H w U' 1 - ihU D d t U ] D . (S16) 

By direct substitution of A = A(t ) + D into the Hamiltonian we find: 

A'A = {A* + D'){A + D) = |H| 2 + (AD* + A*D) + D*D (S17) 

and therefore: 

l h L 

Hw,d = hXo p| 2 + (AD^ + A*D) + Z?l£)] + XjBjBj + — Hi P q S BjB^ p B q B s + he(t)M 00 (D^ + D + A* + A) 

j =1 l P qs=l 

L 

+ 2 h p| 2 + (AD* + A*D) + D*D] VipBjBp - i{A(t)A I - A(t)*A) + he(t) ^ M W {B\ + B t ). (S18) 

Ip 1 = 1 


We choose A(t) by requiring the terms linear in D and D 1 to vanish 

hXoA(t) + (2h'^2 r li P Bj B p )A(t) + he(t)M 00 = —ihA(t), 

ip 


(S19) 
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so that we can eliminate the terms involving the drive. Eliminating the drive is equivalent to moving to a frame where 
the cavity evolves according to the usu al H eisenberg equation A = j^[H,A] — ^MqoA, and therefore solving for the 
classical part A(t) is equivalent to Eq. (S19). Since r]i p <C Aq we neglect that term and consider only: 


XoA(t) + e(t)Moo = —iA(t). 

In the stationary case A(t) = 0 the solution is: 

—e(t)M 0 o 


A(t) = 


mo 


(S20) 


(S21) 


where we have included the correction due to the cavity damping rate k. 

The final Hamiltonian then becomes: 

l l L 

H w , d = HX o [\A\ 2 + D'D] + h X A ' 11 < 2 X fJ-lpqs B}B$B q B a + he{t) J2 + B t ) 

3 — 1 lpqs= 1 

+ 2 h [\A\ 2 + (AD^ + A*D) + D*D] X r lipBjB p . 

i P 


i=i 


(S22) 


Ignoring the constant offset term h\o\A\ 2 , we break the Hamiltonian up into the following pieces. We call Hn the 
term involving only the dressed cavity operator: 


h d = nx 0 D^D- 

Hb the one involving only the dressed qubit operators: 


Bp, — h ~y ' XjBjBj + ^ f^ipqsB^B^B q B s ; 

j =1 lpqs=l 

with Hp b the interaction between dressed cavity and dressed qubits: 


h D iB = 2 h [\A\ 2 + ( Ad f + A*d) + d'd] ^ vi p b}b p . 

Ip 

Finally, we have Hb, drive, the driving term on the dressed qubit operators: 

L 

H b , drive = he{t)Y / M oi(Bj +B t ). 


(S23) 


(S24) 


(S25) 


(S26) 


i=i 


The final Hamiltonian is the sum of all of these terms: 


H cS = Hb + Hb + Hb,b + Hb, drive- 


(S27) 


Cooling The term Hb,b is the one we are most interested in. Depending on the detuning of the drive, this term 
can act either as a cooling term or as a term that induces a state-dependent shift of the cavity resonance and therefore 
allows one to detect the state of the system from a homodyne measurement. It is important to notice that this term 
preserves the U(l) invariance of the Hamlitonian Hb and therefore conserves the total qubit excitation number N, 
thus allowing for cooling within a given excitation-number manifold. 

Let us focus on the cooling process first: we assume the loss rate of photons from the cavity to be high with respect 
to the dynamics involving Hb,b- Therefore, we neglect the terms D and A*D and we call the remaining operator 
the cooling term: 


Kool = 2hA(t)rfJ2 r )ip B i B P = 2hA{t)rf Os¬ 
ip 


(S28) 


Within our approximation, the dissipation makes Vj , ineffective. In Eq. (S28) we have defined the operator: 


Ob = VipBjB p . 

i P 


(S29) 
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If the pump is red-detuned from the cavity, then the cooling operator Wool scatters one qubit excitation to a lower 
energy state by creating a photon with higher frequency than the incoming one (strictly speaking, it scatters a pump 
photon up to the cavity frequency). Therefore the total number of qubit excitations is conserved, but some excess 
energy has been transferred from the qubit array to a cavity photon. 

To estimate the rate of cooling between two states we use the Fermi’s Golden Rule: 


2ir 


r cool - ^ |(^'/|bcool|tE'i)| Ef). 


(S30) 


In the Hilbert space of the dressed excitations (qubits) (g> |cavity), the initial state is |’F i ) = \4>i,0). We consider 
the qubits initially in an excited state ipi with energy E, and no photons in the cavity. The final state instead is 
I'S f) = |where the qubits are in a lower energy state ipf with energy Ef, and the excess energy Ei — Ef is 
carried away by a new photon. Therefore we have: 


r 


O 7 r 

= - -(2hA{t)f Q\0^ B D\ipf, 1)(V>/, 1| D'OsWi, 0 )5{E i -Ef- 




= (2nh)(2A(t)) 2 (4> i \0^f)(^f\0 B m ^<01^11X11^10 )m -Ef- e q ) 


= {2A{t)f\M lf \ 


= (2A(t)Y\M lf \ 


f‘+oo 
! —oo 

/*+ OO 

> —oo 
/*+oo 


Q 


= (2A{t)) 2 \M lf \ 2 I dte^-^iOlD^D^O) 

J —oo 

= {^A(t)) 2 \Mif \ 2 Sdd{^>i - u>f), 
where Mif = (^/|C ) s|V’i) and 

SdD= (w - A c ) 2 + (ft/2) 2 


(S31) 


(S32) 


is the spectral density of the cavity field fluctuations. The square of the classical drive amplitude A(t) 2 is equal to 
the average photon number in the cavity: A(t) 2 = n, so we see that the cooling rate is linear in the photon number: 


Ti-c/ = M M if\ 


(w* - Uf - A c ) 2 + («/2) 2 ‘ 


(S33) 


In the stationary case (Eq. (S21|) A(t) ~ e(f), and therefore the drive power is directly proportional to the photon 
number P ln oc n, giving the result contained in the main text (Eq.(4)) that: T coo i oc P ln . 

The transition rates just derived are perturbative, since the derivation relied on Fermi’s Golden rule. Thus they 
are valid only in the limit that the rate of transitions is small with respect to the decay rate of the cavity: 


11/001 ^ ft* 


(S34) 


Beyond this regime the perturbative description of the cooling process is not appropriate anymore because the effect 
of the V^ ol term will no longer be negligible, and will induce other processes besides the cooling [S3] . 


Stark shift For low intracavity photon number, the frequencies w,; and uif used above in Eqs. (S31) and (S33) are 
simply the bare eigenvalues of the Hamiltonian in Eq. (S24) that we diagonalized to find the eigenmodes. For higher 


photon number, the first term 2h \A\ 2 J2i P VipEjB p of the Hamiltonian H b ,b (Eq. (S251) starts to be significant; this 
term contributes a Stark shift to the energies of the eigenstates. We consider this Stark shift perturbatively and define 
a photon-number dependent frequency as: 


u}i(n) = + 2hn(ipi\0 B \ipi)- 


(S35) 


For the low photon number used in the experiment (at most five photons at steady state in the cavity), this perturbative 
correction is a good approximation to the exact solution. 
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Measurement By performing a homodyne measurement, we can measure the shift of the cavity frequency de¬ 
pending on the state of the array: 


Ao + 2 r]ipBj B p 

ip 


D^D. 


Therefore we introduce the operator \ describing the cavity pull: 

X = 2 E r Hp B i B p = 20 B ■ 

i P 


(S36) 


(S37) 


The expectation value of this operator on a generic state of the array S constitutes the observable that identifies that 
state: 


Xs = 2 (0 B ) S = 2(]T mpB\B p ) s . (S38) 

i P 

We point out an important difference in the behavior of the operator Ob during the cooling process versus during the 
measurement. During the cooling process the pump drive is centered at a red detuning A c = w* — uif given by the 
energy difference between the two states, initial and final, constituting the cooling transition. The operator Ob then 
will induce transitions between states whose energy difference lies within a bandwidth ~ k centered on A c . During 
the readout, instead, the drive is at the cavity frequency, i.e. at zero detuning, so to a good approximation we can 
neglect terms rotating faster than the linewidth n of the cavity. For small lattice sizes and big coupling J, as it is 
the case in our experiment with L = 3, the excitation energies are big compared to n and therefore the operator \ 
will not induce transitions between different states during the readout, and the measurement can then be considered 
quantum nondemolition (QND). 


Long array limit 

For long arrays, L 1, and small hopping strength J, the energy levels will become closely spaced and therefore 
eventually the measurement will not be QND anymore. Once there exist pairs of states with energy difference A E less 
than ~ k, the measurement drive will (generically) cause transitions between them. This same physics will adversely 
affect the cooling process that we have discussed. The cooling process implements the removal of an amount of energy 
A E = H(uji — uif) (u>i > uif) that is well-defined up to the cavity linewidth n. A E k is easily achievable in small 
arrays. In the limit of long arrays, eventually A E ~ k or smaller and both Stokes and anti-Stokes processes become 
possible and the (non-equilibrium) quantum bath no longer has zero effective temperature jS39l . 

Consider the simplest case in which all of the quits have uniform frequency and are coupled only to nearest- 
neighbors. Neglecting initially the anharmonic terms, i.e. imagining a chain of coupled linear harmonic oscillators 
gives the simple Hamiltonian 


L L -1 

#array = E b]bj + HJ + b]b j+1 ). (S39) 

0 =1 J=1 

This chain of oscillators exhibits normal modes with excitation amplitudes varying sinusoidally along the array. Or, 
in the language of excitations as bosonic particles, this is a tight-binding model of free bosons on a lattice, which is 
diagonalized by introducing modes of the form: 


B 


n 



sin {k n j)bj, 


k n = 


L+l 


n = 1, 2,..., L. 


(S40) 


In this form the Hamiltonian becomes 


L 

B'rUYHy k, E [w 0 + 2 J COS {kn^B^Bn. 

n =1 


(S41) 
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Note that unlike in the previous derivation we did not include the coupling to the cavity when diagonalizing the 
quadratic Hamiltonian. Including these terms we now obtain 


H = + 2 J cos(k n )]Bl l B n + hhj c a) a + fi E 9j (M 1 + h )°) 

n =1 j=l 

L L 

— cos (/C 77 ) ] . 0 ^-Z^, H - frw c a f a + H ^ ^ £tti H - B^clJ 


71 = 1 


m— 1 


where 


£m = E * sin (*W')- 

j=i 


The quartic interaction term becomes in this basis 


(542) 

(543) 


(S44) 


where 


Hint = 2 E n EW 
1=1 

= E 

m,n,p,q= 1 


mnpq Brn B n Bp Bq , 


1 —rnnpq 


2 

(i + iF 


L 

E Uj sin {k m j) sin (k n j) sin (k p j) sin (k q j). 
i=1 


(545) 

(546) 


(S47) 


Proceeding as before we can diagonalize the quadratic part of the Hamiltonian in the presence of a drive on the 
cavity and extract the cooling operator defined in Eq. (S28) and Eq. (S29) and hence determine the cooling matrix 
elements M,f . We will not present those details here but rather simply note certain qualitative features of the result 
that can be seen more easily in the present version of the derivation where we have delayed inclusion of the cavity 
coupling until after diagonalization of the array Hamiltonian. If gj varies slowly and smoothly with position then 
will be large only for small m. It follows that the Raman scattering process does not permit large changes of (quasi-) 
momentum and the matrix r]i p in Eq. (S29) will be non-zero only near the diagonal. If only small momentum changes 
are permitted then only small energy changes are permitted and the cooling will be weak. If on the other hand, the 
set of {gj} has low symmetry then the {£ m }, and hence the transition matrix elements Mjj derived from them, are 
not constrained by symmetry and are generically non-zero. 

We now discuss a protocol to cool towards the lower energy eigenstates of a general L-site array. We illustrate the 
basic idea using the single-excitation manifold as an example. Based on the band structure of the ID tight-binding 
model, we know that the single-excitation manifold energies are spread around the bare qubit frequency ujq with a 
total width of 4 J. The levels will be more densely spaced at the top and bottom of the band (see Fig. SI) due to the 
van Hove singularities in the extremes of the tight-binding band structure. 

For simplicity we assume that all cooling matrix elements are generically non-zero. Suppose we start from a high- 
energy state Ei and we set the drive such that the cavity detuning k < A c < 4 J can cool between transitions within 
the bandwidth. As illustrated in Fig. |S1| a cascade of several cooling steps brings the initial energy down to a final 
value equal to Eq = Eq mod (A c ) < A c , at which point no more cooling steps can occur. To lower the energy 
further, we can slowly decrease the detuning A c towards ~ k. We will be able to cool to the lowest energy state if 
J(j^y ) 2 > n. Otherwise this process will cool down the N-particle subspace to an effective temperature of T e g ~ k 

[EM). 


It may occur that the above procedure fails because some particular matrix element in the cascade vanishes. This 
can be overcome by scanning cooling pump detuning A c down and up multiple times to find a transition that allows 
escape from the trapped state. 

The cooling rate calculations used above for illustrative purposes are straightforward to carry out within the 
manifold of single excitations. Higher excitation manifolds will exhibit more complex and interesting dynamics 
because as energy is removed from the system by Raman processes, boson-boson collisions will further relax the 
particle distribution function. The dynamics will depend importantly on whether or not the underlying model is 
integrable. These considerations are well beyond the scope of the present experimental state-of-the-art, but could well 
become important in future realizations of quantum simulators. 











/K 


L- 



FIG. SI: Illustration of the temperature limit for the single excitation subspace of a degenerate L-site chain. To be able to cool 
to the ground state we require > k. 


Three-qubit array used in the experiment 


We now specialize the above results to the case of L = 3 sites. The Hamiltonian H a in Eq. (S41 is: 


Ho = fiJ2 A J b } b j 

3=1 


hJ 


2 

E 

3=1 


(bj+i b j 


+ bjbj+ 1 ) -+■ hJio(b\bo + 


b\b\) + HA c a^a + h g 3 (bjat + 6 jo). (S48) 

3=1 


where we have also included the next-nearest-neighbor coupling J 13 between the 1st and 3rd sites. This Hamiltonian 
has the matrix form: 


Ho = h 


acting on a four-component vector v = (a, 61 , 62 , 63 ) 


/ A c 

9i 

92 

93 \ 

9i 

Ar 

J 

•h 3 

92 

J 

A 2 

J 

\ 93 

Jl3 

J 

A 3 j 

), so 

that 

H 0 

= vt H 


(S49) 


The matrix N defined in Eqs. (S 6 ) and 


(S7) for the experimental parameters given in Section [| is: 

/ -0.986 -0.073 -0.126 -0.077 \ 
-0.163 0.426 0.690 0.561 

0.014 -0.680 -0.158 0.716 
V -0.013 -0.592 0.695 -0.408 


N = 


(S50) 


The rows of this matrix are the coefficients of the new dressed operators (A, B\, B 2 , B 3 ) expressed in the old basis 
(a, 61 , 6 2 , 63 ). As expected A ss a remains an operator that is mostly cavity-like, similarly the Bj operators are 
mostly linear combination of qubits operators. This is ensured by the condition gj/(uij — u> c ) <C 1 for each j. In the 
experimental setup we have: gi^/{u 3 i^ — ui c ) ~ 0.06 and ^ 2/(^2 — uj c ) ~ 0 . 12 . 

If the couplings g\ and 53 were identical, the model would be symmetric with respect to the exchange of b± -H> 63 , and 
the operators A, B\, B 3 would be even under such transformation, while B 2 would be odd, transforming as B 2 —B 2 . 
In the experimental case, g\ and g 3 are close but not identical and similarly, the qubit frequencies ui 1 and W 3 are 
slightly different. This symmetry is thus broken in the experiment, but only weakly. We see this, for instance, in the 
matrix N that is almost, but not completely, invariant under exchange of 61 ■<->■ 63 . 

The matrix M = TV -1 , that according to Eqs. (S 8 ) and (S9) has as elements the decomposition of the original 
operators (a, 61 , 62 , 63 ) in the new basis (A, B 1 , H 2 , H 3 ), is: 


M = 


( -0.986 -0.164 0.014 -0.013 \ 

-0.073 0.426 -0.680 -0.592 

-0.126 0.690 -0.158 0.694 

\ -0.077 0.561 0.716 -0.410 j 


(S51) 









































9 


Inspecting the elements of this matrix confirms that for this range of parameters our assumptions Moo 
Mqi -C 1 are valid. 


The tensor 77 in Eq. (S13) is: 


—2.422 0.227 -1.241 
77 = 2n x ( 0.227 -1.279 0.338 
-1.241 0.338 -2.445 


MHz. 


1 and 


(S52) 


The fact that the off-diagonal elements in the above matrix involving the operator B 2 ( 7 ? 12 and 7723 ) are relatively 


small is a manifestation of the near-symmetry our system exhibits, described above. The operator Wool in Eq. (S28) 


cannot connect states with different symmetry (in the perfectly symmetric case Bi, B 3 would be even, while B 2 odd 
and the elements 7712 = 7721 and 7723 = 7732 would be equal to zero). 

Single-excitation subspace Let us consider first the eigenstates with only a single excitation in the array. For 


these states the nonlinear term in Hamiltonian Eq. (S241 vanishes and we are left with a simple diagonal Hamiltonian: 

L 

H b = H^XjB^Bj. (S53) 

i=i 


A basis B for the Hilbert space of the one-excitation manifold is the following: 

| 1 > = 11 , 0 , 0 ) = b\\g), 

12) = |0,1,0) =h||G), 

|3) = |0,0,1) = B\\G). 


(554) 

(555) 

(556) 

(557) 

We have introduced here the notation |to, n, l) <x (B\) m (B 2 ) n (B' l ) l \G) as a shorthand way of representing a normalized 
state with m excitations in the lowest eigenmode, n in the second mode and l in the third mode. Recall that the 
three modes do not correspond to the three lattice sites; each operator Bj has in general a nonzero component on 


each site, as we found in Eq. (S7|. In this one-manifold basis, the Hamiltonian in Eq. (S53) is diagonal, and therefore 
the one-excitation modes {| Ei), i £ [1,3]} coincide with this basis: \Ei) = 1 1 ) , |_E 2 ) = 12), \E 3 ) = 1 3 ) , provided that 

\ 3 . Expressed in this basis, the operator O b in 


Eq. (S29) has the matrix form: 


O b = 


Ai 

< Aj 

: < 

mi 

7721 

7731 

7712 

7722 

7732 

Vl 3 

7712 

7733 


The cooling rate in this manifold, for example, from \E 3 ) to \E\) is easily evaluated from Eq. (S31) as: 

Te 3 ^e 1 = (2A(t)) 2 |r7 13 | 2 S' £ i £ i(A3 — Ai). 

The cavity pull, or yyshift, corresponding to each state is calculated via: 

XEi 2 {0 B ) Bi Zria 


(558) 

(559) 

(560) 


Two-excitation subspace Let us consider the manifold consisting of states which have two excitations in the 
array. We define a basis for this manifold similar to that used for the single-excitation subspace, noting that because 


eigenbasis for the Hamiltonian: 


|1) = |2,0,0) = -^1?I 2 |G), 

(S61) 

|2> = |0,2,0) = -^H| 2 |G), 

(S62) 

|3> = |0,0,2) = -^ J B] 2 |G), 

(S63) 

|4) = |l,l,0)=B 1 t Bt|G), 

(S64) 

|5) = 0,1,1) = b\b\\G), 

(S65) 

| 6 ) = |1,0,1) = b\b\\G). 

(S 66 ) 
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The Hamiltonian in Eq. (S241 


H b 


o o 

HlpqsBj B^BqBg 

j =1 lpqs= 1 


(S67) 


expressed in the basis B has the matrix form: 


/ 

2Ai + /rim 

1*2211 

1*3311 

TWl 

x/21*2311 

x/21*1311 

\ 


1*2211 

2X2 + 1*2222 

1*3322 

%/2l*1222 

721*2322 

\/21*1322 



1*3311 

1*3322 

2A3 + /I3333 

721*1233 

x/21*2333 

x/21*1333 



721*1112 

V2 1*2212 

V2 1*3312 

Ai + A2 + 2/L*i212 

2/*2312 

2/*1312 



721*1123 

'J2 1*2223 

V2 1*3323 

2l*1223 

A2 + A3 + 2^*2323 

2/*1323 


V 

72 /ini3 

V2^2213 

721*3313 

2l*1213 

2/*2313 

Al + A3 + 2/l 131 3 

/ 


and the operator Ob in Eq. (|S29[): 


/ 2t7h 

0 

0 

72*7i2 

0 

72*7i3 \ 

0 

2*722 

0 

72*712 

72*723 

0 

0 

0 

2*733 

0 

72*723 

72?7i3 

72*712 

72?7i2 

0 

*7n + *722 

*?13 

*723 

0 

72*723 

72*723 

*713 

*722 + *?33 

*712 

V 72*713 

0 

72?7i3 

*723 

*?12 

*711 + *733 / 


(S68) 


(S69) 


To find the eigenstates of the two-excitation manifold we need to diagonalize the matrix: 

Hb = Hb + 2H\A\~Ob, 


(S70) 


that represents Hb plus the photon-dependent Stark shift 2h\A\' 2 Y^ipVip^jB p . The six eigenstates of the matrix Hb 
correspond to the six E-states introduced in the main text, with corresponding eigenfrequencies ej, j = 1,2,..., 6. 
The cooling rate from | Fj) to | Fj) is evaluated from Eq. (S31) as: 

Tfwfi = ( 2A(t)) 2 {F l \0 B \F ] )\ 2 S DD {e l - ej). (S71) 

The cavity pull in a generic Fj state is: 


Xf, 2{F j \Ob\F 0 ) 

K K 


(S72) 
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TABLE I: Cooling rates Ti_>/ in MHz 


\1pi) 

IVt) 

Experiment 

Theory 

\F6) 

|F5) 

0 

0 


|F4) 

11.6 

16.1 


|F3> 

4.2 

9.8 


|F2> 

0 

0.28 


|F1> 

0 

0 

\F5) 

|F4> 

0.4 

0.86 


|F3> 

0.5 

0.75 


|F2) 

5.8 

12.5 


|F1> 

0 

0.2 

m 

|F3> 

NA 

10.6 


|F2> 

3.1 

4.2 


|F1) 

8.4 

10.2 

|F3> 

|F2) 

0.6 

0.66 


|F1> 

11 

20.6 

|F2) 

|F1> 

2.6 

10.3 

|£3> 

|F2> 

0.54 

0.52 

\E3) 

|F1> 

13 

15.5 

\E2) 

\E1) 

0.54 

1.15 


SUMMARY OF THEORETICAL RESULTS 


Cooling rates 


As shown in Eqs. (S31) and (S59), (S71) we can theoretically predict the cooling rate for a specific transition as a 


function of the average photon number n and the detuning of the cooling drive from the cavity resonance: 


= An\M if \ 


(u>i —ujf — A c ) 2 + (k/2 ) 2 


(S73) 


As noted above, the frequencies uii and w/ in the above equation depend on the photon number via the Stark shift, 
an effect that was observed experimentally (see Fig. 3b of the main text). At resonance, the cooling rate per photon 
number is then simply proportional to the square of the matrix element of the cooling operator on the two states 
Mif = (%/jf\OB\ipi)'- 


pres n 

= ~\ M if\ 2 - (S74) 

In Table |T] we show the cooling rates per photon number for different transitions in the two-excitation subspace and 
single-excitation subspace, comparing the theoretical predictions with the rates measured in the experiment. The 
theoretical rates are always within a factor of two of the measured ones, and consistently always higher than the 
experimental ones. This suggests that the single-mode cavity approximation, or others, in the theoretical model are 
underestimating other decay processes which weaken the effective cooling rate accessible experimentally. 


We now refer back to the simple tight-binding model of free excitations hopping on the lattice (that we have analyzed 
above in Eq. (S39) and following), we expect this to apply well to the single-excitation manifold where the nonlinear 
terms are zero. On a three-site lattice the only allowed Fourier modes are: k\ = 7r/4, k 2 = 7t/2 and k 3 = 37t/4; 
these correspond respectively to the momentum of the states: E 3 , E 2 , E\. Given that J > 0, the state with higher 
momentum has the lowest energy (here E\ has momentum k 3 = 37t/4). In the case where the couplings obey g\ = g 3 , 
they have parity symmetry. Hence cooling can be effective between states E 3 and E\ because they have the same 
spatial parity. Conversely cooling is highly suppressed in the case of E 3 to E 2 or E 2 to E\ because those transitions 
require a change of parity. This is clearly confirmed by the cooling rate experimentally measured between the states 
of the E-manifold (see lower part of Table [I]). The cooling rate T e 3 ^e 1 is more than twenty times bigger than T e 3 ^e 2 
and Te 2 ^e 1 - 
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Calculation of Ti Purcell-limited relaxation time 


We calculate the Purcell relaxation rate as the decay rate of a qubit-array state due to the action of the bare cavity 
operator a; this operator destroys a photon and induces transitions to lower states. In order to estimate the rate for 
such a process to occur, we compute the overlap between the \Ej) and | Fj) states and the bare cavity mode a. From 
Eq. (S7), we know the decomposition of each dressed qubit operator Bj in terms of the bare operators: 


Bj — Njoa + ^ ' Njibi. 


(S75) 


1=1 


Since the single-particle states \Ej) are simply given by the dressed operator Bj acting on the vacuum (with no 
excitations in either qubit or cavity), the overlap between the | Ej) state and the bare cavity mode is simply given by 
Njo . Then, the Purcell-limited decay rate for this process is the square of this overlap times the cavity decay rate k: 

R Ej ^ g = \{G\a\Ej)\ 2 k = Nj 0 n, T x = — (S76) 

tlR,j=rG 

The results for the | Ej) — > |G) decay times are shown in the top three rows of Table [TT) 

The analysis for the decay rates of the | Fj) states is similar with the decay rate of, say, an initial two-excitation 
substate to a single-excitation state given by 


RFj=fEi = \{Ei\a\Fj)\ 2 K 


(S77) 


_ jp . 

To calculate these values we expand | Fj) in the basis B of states {|n)}„ e r 1 ,6] in Eqs. (|S61j)-(IS66J) : | Fj) = E„=i c « J | n )- 
Then the decay rate to a specific one-excitation state | Ej) is 



6 

2 

6 

RFj=,Ei = 

(Ei\aJ2 c n\ n ) 

71=1 

K = 

J2 c n(Ei\ a \n) 

n—1 


(S78) 


and the decay rate for a direct decay to the ground state given by 



6 

2 

6 

Rf^G = 

(G|a 2 ^c n | n ) 

71=1 

K = 

^]c n (G|a 2 |n) 

71=1 


(S79) 


The total Purcell-limited decay rate is the sum of the rates given above, with a Purcell-limited Tf time given by 
the inverse of that sum. The T\ times found from the above analysis are shown in Table [TlJ For comparison, the Tf 
decay times measured in the experiment are shown in Tabic |III| for the downward rates. Due to spurious excitation 
from the sample being at finite temperature, the experiment observed upward transition rates as well. These spurious 
excitation rates are shown in Table EY1 

In most of the cases, the theoretical prediction for Tf is of the same order of magnitude as the experimentally 
measured value. Nevertheless, there are some qualitative discrepancies. For instance, the theoretical Tf times for the 
E-states are in general shorter than the measured ones. We attribute this to the limitation of the single-mode cavity 
model which can be shown to predict a shorter life-time than a model which takes into account higher modes of the 
3D cavity. A more mundane discrepancy arises from the fact that as the \E X ) and |£? 2 ) are almost dark states, their 
dominant decay channel is not Purcell decay via coupling to the cavity, but rather on-chip material losses. 


Dark and bright states 

Coherently preparing an array eigenstate is accomplished via pulses applied to the same port of the cavity used to 
perform cooling and state measurement. Thus, the response of the system to this drive is contained in the operator 
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TABLE II: 7\ theoretical Purcell-limited decay time in gs 


state 

u/(2n) GHz 

Ti(tot) 

Tx(Ex) 

Tx(E 2 ) 

Tx(E 3 ) 

Tx(G) 


4.61164 

97 




97 

e 2 

4.85539 

80 




80 

E 3 

5.0196 

0.6 




0.6 

Fx 

9.11862 

20 

32 

201 

93 

438 

f 2 

9.3201 

7.5 

0 ° 

oo 

57 

705 

>lms 

f 3 

9.48676 

1.3 

1.3 

50 

212 

>lms 

Fx 

9.64465 

1.2 

1.3 

69 

30 

182 

Fs 

9.7987 

0.6 

49 

0.6 

159 

>lms 

f 6 

9.97278 

0.9 

78 

>lms 

1.17 

5.7 


TABLE III: Ti experimental fitted downward decay time in gs 


state 

uj/(2tv) GHz 

Ti (tot) 

Tx(Ex) 

Tx(E 2 ) 

Tx(E 3 ) 

Tx(G) 

Ex 

4.6078 

28.5 




28.5 

e 2 

4.7854 

30.5 




30.5 

e 3 

5.06916 

3.2 




3.2 

Fx 

9.1184 

18.0 

20.4 

151 



F 2 

9.2592 

15.1 

30.3 

30.1 



f 3 

9.4230 

OO 

bo 

13.5 

25.4 



Fx 

9.5618 

4.6 

5.3 

34.6 



F 5 

9.7788 

3.1 

100.7 

3.3 

70.0 


f 6 

10.0539 

1.5 

20.6 

42.6 

1.6 



Hb, drive discussed in the previous section (see Eq. (S26) and its surroundings). In the dressed basis, this operator 
is Hb, drive = fte{t ) Moi (B\ + B{), so the relevant value(s) of M 0 i for a particular eigenstate determine the 
magnitude of the response. For the single-excitation manifold, this is an easy calculation, as the operators Bj create 
a single-particle eigenstate from the vacuum. For the two-particle manifold the calculation is more involved, as the 
eigenbasis of BjB i does not coincide with the two-particle eigenmodes. A second, equivalent method is to work in the 
bare basis and compute directly the matrix element of the operator H lnt = h ■_ 1 g :) (bj (d + b^a) between the states 
|S, 0 p h) and |G, l p h), for desired array eigenstate | S). This calculation yields 


ds.G = | | Flint |'ll g) I = ft 


(0 P h|a|l ph )(5|^ ft 6t|G) 


3 =1 


(S80) 


In Fig. S2 we plot this coupling <Ie,g for the |l£)-states to ground state G, as a function of flux (current in the coil). 
The figure shows that the theory is qualitatively in agreement with the measurements. | E 3 ) is predicted to be always 
bright, in agreement with the spectroscopy image Fig S4 Both \E\) and \E 2 ) states instead become dark at a specific 


value of the flux (current in the coil). There is some uncertainty in the location of this dark spot due to the uncertainty 
A gj = ±7MHz in the measured values of the couplings g :i , which affects this result significantly. Theoretically we find: 
fdark(Ei) = 11.3 ± 0.7nrA and Ida.A(E 2 ) = 13.1 ± 0.7nrA. These predicted values are shown, with the corresponding 
error bar, in Fig. S2 Experimentally the measured values are: ^ark^i) = 10-71 mA and I^ k (E 2 ) = 10.64mA; these 


are marked with a single dashed vertical line in Fig. S2 (the two lines are too close to be distinguished on the scale 


TABLE IV: Ti experimental fitted upward decay time in gs 


state 

co/(2n) GHz 

Ti(tot) 

Tx(Ex) 

Tx(E 2 ) 

Tx(E 3 ) 

Tx(G) 

Ex 

4.6078 

82.4 




82.4 

e 2 

4.7854 

182.3 




182.3 

e 3 

5.06916 

167.0 




167.0 

Fx 

9.1184 

104.2 

104.2 




f 2 

9.2592 

82.0 

237.5 

125.2 



f 3 

9.4230 

45.3 

98.1 

84.2 



Fx 

9.5618 

36.7 

43.3 

239.4 



Fs 

9.7788 

9.7 

50.8 

28.2 

20.7 


Fe 

10.0539 

3.0 

9.8 

33.6 

5.0 
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of the graph). The measured value (£d) fits within the uncertainty range of the theoretical prediction; however, 
^dark(^ 2 ) is somewhat off the theoretical prediction, and it is also slightly smaller than //)//.(E)), contrary to what 
the theory would predict. 



FIG. S2: Bright and dark features of the 15-manifold states. The figure shows the predicted coupling ds,G in Eq. (S801 between 
a state S in the single-excitation manifold and the cavity pulse. The vertical line show the experimental value for which the 
states Ei and E 2 become dark (the values are very close so cannot be distinguished). The error bars for the theoretical location 
of the dark point for the E\ and Eo states are shown along the horizontal axis. 


EXPERIMENTAL DETAILS 

Device Parameters 

Our device consists of three transmon qubits |S33| on a single silicon chip. Each qubit is formed by two aluminum 
paddles, connected by either a single double-angle-evaporated Al/AlO^/Al Josephson junction (middle qubit) or a 
superconducting quantum interference device (SQUID) consisting of two junctions (outer qubits). The qubit array is 
located in the center of a copper waveguide cavity with dressed frequency ujJ 2it = 7.116 GHz and k/2it = 10 MHz. 
Each qubit couples to the cavity via a Jaynes-Cummings a x (a + a' ) interaction with strength gi. Because the middle 
qubit is located in the center of the cavity where the E-field strength is greatest and its paddles are longer, it couples 
to the cavity more strongly than the outer qubits, with strengths g m id/27r = 264 ± 7 MHz and g 0 ut = /2n 155/149 ± 
7. The outer qubits are characterized by a charging energy E c /h = 214 MHz, which also gives us the absolute value 
of the anharmonicity for a transmon (a = —E c /h), and have a Josephson energy which gives, at zero flux, w q i/27r = 
5.074 GHz for the left qubit and oj q z/2t: = 5.165 GHz for the right. For the middle qubit, E c /h = 240 MHz and the 
qubit frequency is w q 2/27r = 4.892 GHz. The qubits are spaced by 1 mm, giving a nearest-neighbor coupling strength 
of J/h = 177 MHz and a next-nearest-neighbor coupling stength of = 26 MHz, with uncertainties of a couple of 
MHz mainly due to the uncertainty in the calibrated values. The qubit-cavity couplings gi and the qubit charging 
energy E c /h were independently calibrated by suppressing the qubit-qubit interactions. To do so we attached two 
coils to the cavity. One, wraped around the whole cavity, gave a uniform field. The other was fixed off center on top 
of the cavity, which produced a gradient field. Using the combination of these two coils we tuned the qubits such 
that they are effectively not interacting, A q . qj > 10 J qiqj - The qubit frequencies at zero flux as well as the qubit-qubit 
couplings were obtained by fitting the spectroscopically-measured eigenergies of the ten lowest-lying eigenstates of the 
array to the Bose-Hubbard Hamiltonian, as described in the next section. 
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^0.29mrn 

FIG. S3: To scale layout and dimensions of the chip with the three transmons. Josephson juctions are not illustrated 


The dimensions and layout of the chip placed in the cavity are shown in Fig. S3 As can be seen, the transmons 


have slightly different dimensions, so that they interact with the cavity with different strengths (since the field of 
the fundamental cavity mode is roughly uniform over the dimensions of the small chip, the difference in interaction 
strengths comes primarily from the different antennae configurations). In our chip the interaction between the cavity 
and the transmon in the middle is nearly twice as strong as that between the cavity and the transmons on the end of 
the array (recall from the spectroscopy results presented in the main text that the measured <72 was 264 MHz while 
the measured g ±3 was 150 MHz). Numerical simulations indicated that this mismatch in coupling strengths would 
improve achievable cooling rates. 

Coupling between the transmons themselves comes from two sources: cavity-mediated interactions and direct 
capacitve (or dipole-dipole) coupling. Cavity-mediated interactions arise when both qubits couple to the cavity mode, 
as described in ref. This coupling strength is Jj 7 cav it y = \gi9j{ 5 - + ^-), about 20 MHz for adjacent qubits 

and 10 MHz for edge-edge coupling, at the working point of the experiment. Direct dipole-dipole coupling, discussed 
in refs. [S36 1 IS401 . arises from the capacitance between transmon paddles. In our case, this coupling is on the order 
of 150 MHz for adjacent qubits, and 10 MHz for the qubits on the edge. 


Spectroscopy 


We perforin spectroscopy to extract, as a function of applied magnetic flux, the eigenergies of the nine lowest-lying 
excited states of the array with respect to the global zero-particle ground state. These nine states consist of three 
states in the single-particle manifold and six in the two-particle manifold. As in the main text, these states are denoted 
|G), {| Ei}}, and {|F))} respectively. We probe the array for coil currents between -2 and +17 mA. 

Since the qubit population without any excitation predominantly lies in |G), standard two-tone spectroscopy reveals 
the |G) -+ \Ei), |G) —>■ \E 2 ), and |G) —► \E 3 ) transitions. To perforin this spectroscopy, the reflected phase of a tone 
near the cavity resonance (7.116 GHz) is continuously monitored as a second tone sweeps from 3.7 to 5.3 GHz. This 
measurement results in Fig. |S4[ i, with three main lines indicating the single-particle energies. 

Extraction of the two-particle energies is more involved. For the \Fq) state, the energy can be directly measured 


via a two-photon transition from |G), as shown in Fig. S4i. For all other |Fj) states, however, the energies must be 


measured indirectly via transitions from a single-particle state. We use | E 3 ) and \E 3 ) as stepping stones to measure 
the | E 3 ) —► | Fi) and \E 3 ) -+ Iq) transitions, by running two additional spectroscopy scans: one with the addition 
of a tone at the |G) -+ | E\) frequency, and another with the addition of a tone at the |G) —> | E 3 ) frequency. The 
results, shown in Fig. |S4 ) d, allow the identification of all six two-particle states. 

From the extracted energies of the one- and two-particle manifolds (see Fig. 1 in the main text), we extract 
parameters of our device by fitting these values to predictions based on the Bose-Hubbard Hamiltonian with an 
additional next-nearest-neighbor coupling term. After taking into account the variation of the qubit frequencies with 
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Coil Current (mA) Coil Current (mA) 


FIG. S4: Raw data from spectroscopy, showing (a) the |G) —► |£)) transitions probed with two microwave tones, and (b) the 
l-Ei) —> | Fi) (blue) and I-E 3 } —> \ Fi) (red) transitions probed with three microwave tones. 


flux, that Hamiltonian is 

3 

H = hJ2 (ui(<f>)b\h + yStft&Si) + hJ (b\b 2 + b\b 3 + Ac.) + J 13 (b\b 3 + Ac.) (S81) 

2=1 

with the qubit frequency as a function of flux described by w(</>) = woii/cos (13,1 + A) for the outer qubits with 
SQUIDS (for the middle qubit u is constant). The parameters Bi for each edge qubit are the ratio between the 
current applied to the coil and the flux threading the qubit’s SQUID loop, and A is an overall offset due to potential 
flux trapped in the SQUID loops during cooldown. 


Dispersive readout of the array 


Due to the dispersive coupling between the qubits and the cavity, each array eigenstate induces a shift in the resonant 
frequency of the readout cavity. The frequency corresponding to the array in |G) is measurable simply via the reflected 
phase measurement on a network analyzer, as the array is in its global ground state without excitation pulses. This 
frequency is 7.116 GHz when our system is biased up to 10 mA. To measure the resonator frequencies corresponding 
to the excited states, we use microwave pulses to prepare the array in the desired eigenstate |i), then measure the 
reflected phase 9i of a 7.116 GHz tone, referenced to the reflected phase with the array in |G). This measurement was 
done using our LJPA in phase-preserving mode. The standard equations for a reflected phase shift from a resonator 
yield that the frequency shift Xi for a given eigenstate is related to 0* by the equation Xi = k/ 2 tan (0,/2). The 
measured reflected phase angle 0 exp and the corresponding y exp are shown in Table [V] In the same table are also 
shown the x shifts calculated theoretically according to Eqs. (S60l and (S72). Except for \F 3 ) and Aj), all of the 
states are resolvable. In fact, by using the LJPA in phase-sensitive mode and adjusting the measurement frequency 
and amplification axis (phase of the detected quadrature), we can obtain additional separation between states of 
interest for a particular experimental run. 

To extract the population of a given state after a particular experimental sequence, we repeat the sequence several 
(~ 1 million) times and histogram the measured phase or quadrature amplitude values. After the run, we take a set 
of calibration histograms, in which we prepare all ten states and immediately make a measurement with the same 
frequency and amplification axis used in the experiment. We then fit the measured histograms to a sum of Gaussians 
with the same mean and variances as the calibration histograms, and from the amplitudes of each Gaussian, extract 
the corresponding state’s population during that run. See the next section for an example of this procedure. 
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TABLE V: x shifts 


state 

#ex P (rad) 

Xexp/ K 

Xth/« 

El 

1.37 

0.41 

0.49 

E2 

0.74 

0.19 

0.26 

E3 

1.43 

0.43 

0.48 

FI 

2.09 

0.86 

1.07 

F2 

1.64 

0.53 

0.68 

F3 

1.82 

0.64 

0.75 

F4 

1.77 

0.61 

0.70 

F5 

2.03 

0.80 

0.82 

F6 

2.16 

0.93 

0.90 



FIG. S5: An example of calibration histograms. 


Extraction of Decay Rates 

Natural decay dynamics were measured by coherently initializing the system in the state of interest, waiting for a 
time, then measuring the array state, as described above. We fit the population dynamics for each state to a decay 
model given by the matrix rate equation <9 t c = Tc, where c is a vector containing all of the state populations as 
a function of time and T a matrix with transition rates between states. This procedure is similar to that used by 
Peterer et al. |S38I in their study of a single transmon qubit. This model suppresses transitions from the two-particle 
manifold directly to the zero-particle ground state based on the results in Ref. |S38j . We also suppress transitions 
between states in the same manifold, since these transition frequencies are on the order of a few hundred MHz, and 
the photon shot noise spectrum which plays a dominant role in dissipation for this system has very little support at 
these frequencies. Best-fit parameters are given in m and were used to generate the natural decay map in Fig. 2 of 
the main text. For the bath engineering decays, a similar model was used for the fit, with the addition of parameters 
for the intramanifold decay; all intermanifold rates were held fixed to the previously-measured natural values. 

Errors in the fit, primarily at low population, occur likely due to effects such as spontaneous T\ decay which cause 
the readout histograms to be skewed from their nominally Gaussian shape. See |S42I for a detailed explanation of 
this effect. As can be seen in the main text, this effect seems to occur more often in the case of the natural decays 
than during the cooling protocols. 
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Calibration of Bath Engineering Drive Power 


Calibration of powers used to drive the bath-engineering transitions is done in the standard circuit-QED man¬ 
ner jS43| : first the qubit-cavity ;y-shift is calibrated, then with this value known, the intracavity photon number 
is then inferred from the Stark shift on the qubit. To do this, in our system, we tune the edge qubits to below 3 
GHz, decoupling from the middle qubit. We then measure the x-shift between the middle qubit and the cavity by 
measuring, for a variety of incident powers at the cavity frequency, both the measurement-induced dephasing rate 
and the resulting Stark shift. As shown in Ref. [S43] . the measurement-induced dephasing rate is 8 x 2 n/n, while the 
Stark shift is 2 \n, so by comparing the slopes of these lines vs. power, we extract x (since k is known). From there, 
we use the Stark shift to calibrate the intracavity photon number over the range of frequencies and powers used in 
the bath-engineering experiment. 


Comparison of Cooling and Coherent Driving to Excite a Dark State 


Here we include a brief description of our attempts to populate dark states using coherent microwave drives rather 
than via the autonomous feedback protocol discussed at the end of the main paper. 

In theory, without dephasing and dissipation, and in the absence of higher Jaynes-Cummings energy states, one 
could drive coherent transitions to the desired state except at the singular flux bias point where the matrix element of 
this transition is exactly zero. In practice, for experimentally available drive powers, driving an almost-dark transition 
may take such a long time that dissipation (qubit T), T^, cavity loss) reduces the fidelity to unacceptable values. This 
was the case at a flux bias of 10 mA (where the single-particle states were almost dark), where, in contrast to the 
usual coherent pulses which take tens of nanoseconds (we used 64 ns for the coherent pulses), the pulse took around 
1 fjs to excite the |G) —> \ E{) transition, and this caused our fidelity to be limited to about 65% (compare the cooling 
process, which achieved a fidelity for the single-excitation state of 80%). 

Another limitation of the coherent excitations is that, due to higher levels in the spectrum, off-resonant transitions 
which have a higher dipole moment than the dark transition may be driven before the desired transition. This was 
the case at 10.71 mA (“complete” darkness), where even at the maximum power we could excite with, we saw no 
population in the \Ei) state, but at these powers higher states (either in the two-excitation subspace or in an even 
higher manifold; from the dispersive shifts it was difficult to tell properly) began to get populated. 
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